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C"^ Abstract 

A finite difference sclieme is presented for tlie Dirac equation in (1+1)D. It can liandle space- and time- 
rs dependent mass and potential terms and utilizes exact discrete transparent boundary conditions (DTBCs). 

Based on a space- and time-staggered leap-frog scheme it avoids fermion doubling and preserves the disper- 
sion relation of the continuum problem for mass zero (Weyl equation) exactly. 

Considering boundary regions, each with a constant mass and potential term, the associated DTBCs are 
derived by first applying this finite difference scheme and then using the Z-transform in the discrete time 
variable. The resulting constant coefficient difference equation in space can be solved exactly on each of the 
two semi-infinite exterior domains. Admitting only solutions in I2 which vanish at infinity is equivalent to 
imposing outgoing boundary conditions. An inverse Z-transformation leads to exact DTBCs in form of a 
I convolution in discrete time which suppress spurious reflections at the boundaries and enforce stability of 
^ the whole space-time scheme. 

Q An exactly preserved functional for the norm of the Dirac spinor on the staggered grid is presented. Sim- 

O ulations of Gaussian wave packets, leaving the computational domain without reflection, demonstrate the 

!/3 quality of the DTBCs numerically. 

O 

' ^ Keywords: Dirac equation, flnite difference method, leap-frog scheme, dispersion preservation, 

transparent boundary conditions 
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T— I 1. Introduction 
> 

t We start out with a brief summary regarding important properties of the Dirac equation, its role played 

00 in physics, existing numerical schemes for .ts solution, and the issue of open boundaries. 

1.1. The Dirac equation 

Next to its fundamental role in relativistic quantum mechanics and field theory, which provide the foun- 
^y-^ dation of modern nuclear and high energy physics [1] [5] , the Dirac equation has received a rapidly growing 

importance in condensed matter systems as well. Especially, in the context of the recent experimental real- 
L| ization of graphene [3] , 2D and 3D topological insulators [H E] , and optical lattices [BHS| the Dirac equation 

. ^ describes the underlying physics as an effective field theory. Historically it was proposed by Dirac with his 

ingenious idea of linearizing the square root of the relativistic energy momentum relation by the introduction 
^ of Dirac matrices and multi-component wave functions, nowadays known as Dirac spinors. Imposing the 

condition that the twofold application of this Dirac operator onto the spinor must yield the Klein-Gordon 
equation, leads to the Clifford algebra for the Dirac matrices. In (l-l-l)D and (2-|-l)D the minimum dimen- 
sion for a representation of this group is two, whereas in (3-|-l)D the Dirac spinor must have a minimum 
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of four components. However, one can work with two components if one accepts higher-order derivatives 
in space and time [U [2] . The latter has lead to the prediction of anti-matter and the concept of the filled 
Fermi sea in many-particle physics. Alternatively, the Dirac equation emerges from an investigation of the 
transformation properties of spinors under the Lorentz group [S]. In condensed matter physics Dirac-like 
equations arise in context of low-energy two-band effective models, e.g. in k ■ p perturbation theory or the 
tight-binding approximations [5] . Indeed, the study of Dirac fermion realizations has developed into one of 
the most exciting current topics of condensed matter physics [S]. 

In this work we restrict ourselves to a study of the (l-l-l)D Dirac equation and present a numerically stable 
scheme for its solution under open boundary conditions. The latter are motivated by a particle transport 
situation, as well as the fact that, unlike for the Schrodinger particle, a deep one-particle potential does 
not ensure confinement for Dirac particles. In its Schrodinger form, also called the standard or Pauli-Dirac 
form, the (1-|-1)D Dirac equation may be written as (using c = 1, e = 1, and h = 1) 

idt'tp{x,t) = Hip{x,t) , H = m{x,t)a2 — idxCTx — V{x,t)l2 ■ (1) 

The (Ti's are the 2x2 Hermitian anti-commuting Pauli-matrices, and I2 is the identity matrix, x G R, t G R~^, 
and i/" = (u, v) € is the complex 2-spinor. m, F € IR represent, respectively, a space- and time-dependent 
mass and scalar potential. For constant coefficients, Fourier-transformation in the space and time variable 
and solving the eigenvalue problem gives the energy spectrum E± = ±-y/ + p^. In the physical ground 
state all negative energy states are filled with fermions. Empty negative energy states are reinterpreted as 
filled hole states (anti-particles) with positive energy, in analogy to multi-band systems in (non-relativistic) 
condensed matter physics [T]. The norm of the spinor is defined as \ip{x,t)\2 = A/|M(av^OP^R^K^^^^)l^j and 
its square can be interpreted as the probability density in space for given time t. The norm is conserved 
because Eq. ([T]) is of Schrodinger form and H is Hermitian [T]. Global gauge-invariance holds as for the 
case of the Schrodinger equation: the addition of a constant V simply adds a phase to the solution. The 
free Dirac equation has various other symmetries [I] [10] . First there are the continuous transformations of 
spatial rotation (only meaningful for more than one space dimension) which, with the Lorentz boost, form 
the Lorentz group. Together with space and time translations, they form the Poincare group. The latter 
is the fundamental group in particle physics, but is of minor importance in solid state physics because the 
crystal lattice necessarily breaks these continuous symmetries. The discrete symmetries holding for arbitrary 
constants me R,V G R are space refiection (parity) and time reversal symmetry which will also be present 
in the finite difference scheme. 

1.2. Numerical aspects 

Several schemes have been proposed and used for particle transport simulations based on the time- 
dependent Dirac equation. For a numerical treatment one has to discretize the continuum problem either in 
real- or Fourier-space, or a combination thereof. Real-space schemes are, for example, the finite-difference 
[TTVfT^ and finite-element methods |T3], whereas the spectral methods [T5] are examples for the momentum 
space approach. Split-operator methods separate the time-evolution operator into several parts, with each 
of them depending only on either momentum or position |161 117j . While having the advantage of a natural 
implementation of space- and time dependent potential and mass terms, the finite-difference and finite- 
element schemes have to deal with the issue of fermion doubling, which means that, for a given sign of 
the energy, there are two (or more) extrema in the m 7^ energy-momentum dispersion relation instead 
of the single one of the continuum problem |12j . In fact the 'Nielsen-Ninomiya no-go theorem' forbids the 
existence of a single fermion flavor for chirally invariant fermions on a regular grid without breaking either 
translational invariance, locality, or Hermiticity [TS]. In (1+1)D one can get rid of the fermion doubling 
using a staggered grid for the two spinor components. This is equivalent to taking the left-sided first-order 
derivative operator for one component of the spinor and the right-sided for the other one |12| . One obtains 
a monotonic dispersion relation with only one minimum. Here we present a scheme which provides an even 
better result by applying staggering to both the space and the time coordinate. This yields a numerical 
scheme which preserves the exact dispersion relation of the continuum problem for the special choice of the 
ratio between time and space grid r := At/ Ax = 1 and mass m — 0. For rrijV ^ and r = 1 the dispersion 
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relation improves for all possible wave-numbers k e [— tt/Ax, Tr/Aa;] with the refinement of the grid. For 
Ax the numerical dispersion relation becomes identical to the continuum one. This is not true for most 
finite difference schemes in general and, to our knowledge, no finite difference scheme with this property for 
the Dirac equation has been reported before. 

1.3. Boundary conditions 

For a numerical treatment of a differential equation, such as the Dirac equation, the number of degrees of 
freedom must be finite. In addition to the discretization of the time and space variable in a real-space scheme 
one has to restrict the simulation domain to a finite region in time and space. Then, appropriate boundary 
conditions are needed to ensure that the solution obtained within the finite domain is (at least) a good 
approximation to the solution of the whole space problem. Generally, the time-dependent Dirac equation is 
solved as an initial-value problem in time. The standard approach for the derivation of spatial transparent 
boundary conditions (TBCs), e.g. for the Schrodinger equation, has been to solve the continuous exterior 
problem by using the Laplace-transformation in time and to discretize the continuous TBCs afterwards [TO] . 
To avoid stability problems j20| and spurious reflections at the boundary from inconsistent discretization 
schemes, recently, a new improved approach in which the whole domain is discretized first and then solved 
exactly in the constant-coefficient exterior domains using a Z-transform in time has been developed for the 
Schrodinger equation [3TJ [21] . In this way one maintains the stability properties of the scheme on the whole 
space and avoids inconsistent discretization for the simulation region and the boundary conditions. The 
resulting discrete transparent boundary conditions (DTBCs) which are non-local in time are exact in the 
sense that they do not introduce a procedural error. An alternative approach which, in most cases is easier 
to handle because an inverse Z-transform can be avoided, is to discretize in space only, derive TBCs, and 
discretize in time thereafter. This has been done with good results for hyperbolic systems [25H27] . However 
for the proposed scheme, where the good properties regarding dispersion and conservation of norm arise 
from the simultaneous discretization of space and time in an interlaced manner, paying the extra price in 
form of an inverse Z-transformation is well justified. Therefore we follow the fully discrete approach of [21] 
to develop TBCs for the time- and space-staggered leap-frog scheme. Furthermore, in order to preserve the 
covariant symmetry of the Dirac equation on a space-time grid, time and space coordinates must be treated 
on an equal footing. This is particularly important in higher dimensions. 

As a major difference to the Schrodinger equation one should mention that a massless relativistic particle 
cannot be trapped by a scalar one-particle potential [1]. The same holds for massive relativistic particles 
with energy E, when the potential depth V is such that E + mc^ > V > E — mc^. This phenomenon 
is related to the Klein paradox and is due to the two-band-nature of the dispersion relation, consisting 
of an electron and a positron band. One is necessarily dealing with a scattering problem whenever there 
are no bound states supported by the Hamiltonian. Then, for initial data which are compactly supported 
on the computational domain, the solution always reaches the boundary after a finite time. For simulation 
times beyond that threshold, open boundary conditions are required to close the finite difference scheme in a 
particle transport situation, similar to non-relativistic particle transport simulations in nano devices \22\ 123] . 



2. Continuous transparent boundary conditions for the Dirac equation in (1+1)D 

In this section we derive continuous TBCs for the Dirac equation Eq. ([T]). We divide the entire space 
into the computational domain (0, L), and the semi-infinite exterior domains {—oo, 0] and [L, oo). The mass 
is assumed to be constant m{x,t) — m and the potential is constant in space V{x,t) — V{t) in the exterior 
domains. This reduces to the case V = with the following gauge change of the spinor: 




(2) 
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With X — ('^i the exterior, we have 

u{x^ i) 



u{x, t) 
v{x, t) 



v{x,t) ^ ^ 

The multiphcation with ax and a Laplace transformation with respect to the time variable leads to 



u{x, s) 
v{x, s) 



—s + im 

-s — im 



Eq. Q has the general solution 



u{x, s) 
v{x, s) 



-s + im 



— + m"^ 



C2 



ii{x, s) 
v{x, s) 



(3) 



(4) 



(5) 



where is written for the square root with positive real part. Because the solution must be in L2{R) 
the constant C2 must vanish on the right exterior domain. For the same reason ci must be zero on the left 
exterior domain. Therefore, the boundary conditions on the right boundary are 



dxu{x, s) |^_^ = — + m'^ u{L, s) and dxv{x, s)\^_j^ — — \/s'^ + v{L, s) . (6) 

On the left boundary one gets 

dxu{x, s)\^^^ = X/s"^ + t2(0, s) and dxv{x, s)\^^^ = X/s"^ + m^ v{0, s) . (7) 

The structure for both spinor components and both boundaries is the same, so we proceed with u on the 
right boundary. First we derive boundary conditions as a Neumann-to-Dirichlet map by writing 



s) 



1 



-dxu{x,s)\ 



x—L 



Then the inverse Laplace transformation leads to the convolution at the right boundary: 
t) = -Jo{mt) *t dxu{x, l^^i = - / Jo{mT)dxu{L, t - t) dr , 

JQ 

with Jo being the Bessel function of first kind. Analogously for the left boundary: 

u{0,t) = Jo(mt) *t dxu{x,t)\^^^ . 



(8) 



(9) 



(10) 



Finally, this leads with the gauge Eq. ^ for y 7^ to the TBCs in the form of a Neumann-to-Dirichlet 



map: 



iP{L,t) = -e^^-«{jo(mi) n [dxip{x,t)\ 

^{0,t) = e^^'W{jo(mi) *t [dxiP{x,t)\^^^e-'^^'^'^] } 



I ... right TBC , 
. . left TBC . 



(11) 
(12) 



For the derivation of TBCs in the form of a Dirichlet-to-Neumann map we write: 

dxu{x,s)\^^^ = -9{s)u{L,s) , 

with: 

/ / s"^ -\- Tfi^ \ 

g{s) := V + m2 = - s + s . 

\ V S + Tt^ / 



(13) 
(14) 



We use 



and 



fit) = Jo(mt) 



This leads to the following inverse Laplace transform of g{s): 
~g{s) ^-^ g{t) ^ {4{mt) + Jo(mi)] + 8' 



{J2{jnt) + Jo (mi)] +(5' 



(15) 
(16) 

(17) 



where 5' is the first derivative of the delta distribution. Then the TBC as a Dirichlet-to-Neumann map on 
the right boundary is 



On the left boundary one gets 



(18) 



\^\.h(mt) + Umt)\ V(0,Oe-'^'-('^ +att^(0,Oe-*^'^')} . (19) 



Discretizations of Eqs. (11) and (12) or Eqs. (18) and ( |19[ ) can serve as boundary conditions for arbitrary 
finite difference discretizations of the Dirac equation Eq. ([ij. But one has to be aware of the fact that 
inconsistent discretization of the differential equation and the associated boundary conditions usually leads 
to spurious reflections or even instability. As already mentioned in the introduction we will therefore first 
apply the discretization scheme to Eq. ( II for the the boundary regions also and derive the associated TBCs 
by Z-transformation. Eqs. (11), (12), (18), and (19) serve as a guide to gain intuition for the behavior of 
the convolution coefficients. 



3. Time- and space-staggered leap-frog scheme 

Leap-frog time-stepping, in combination with a staggered spatial grid, plays a special role among the 
finite difference methods because, in addition to the elimination of the fermion-doubling problem, it proves 
to be dispersion relation preserving in (l-l-l)D for the 'golden ratio' of r = At/Ax = 1, m = 0, and V — ^ 
(Weyl equation). Moreover, for to 7^ and/or F 7^ the dispersion relation is still monotone and improves 
with a refinement of the grid. It is identical to the exact analytic dispersion when At, Ax — ^ 0, for fixed 
r = 1. This is particularly important for simulations where the whole possible range of the wave numbers 
^ ^ [~Sr' Se] used, for example, in strong external fields. An initial wave-packet which consists only 
of wave-components near fc = may acquire high wave number components due to strong spatial and/or 
temporal changes of potential and/or mass. 

For schemes with non-monotonic dispersion a problem can also arise at the boundary. The modes of such 
a scheme consist of additional, spurious, numerically generated modes on the lattice which, nevertheless, 
must be accounted for in the DTBCs for consistency. This requires special attention at the boundary because 
improper realization of the boundary condition may lead to 'energy' transfer between the modes, spurious 
reflections, and eventually instability (see e.g. [17]). A correct dispersion relation means correct phase and 
group velocity on the grid and is essential for faithful long-time propagation studies. We now present such 
a scheme. 
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3.1. The discretization scheme 

We shall consider the following leap-frog discretization of the Dirac equation in Pauli-Dirac form given 
in Eq.([l]): 



n+l/2 n-1/2 n+1/2 , n-1/2 

U- — U- U- +U 



^ +^im--Vp)^ ^ + ^^-0' (20) 

^n%-^^l/2 ^n+1/2 , ..»+l/2. ^n%+-"-l/2 , (W-+V^),-_V2 _ „ 

At H"^j-i/2 + ) 2 + Ax ~ ' ^ ' 

with j S Z, n € Kl. Here we used the notation rp{x,t) = [u{x,t),v{x,t)] with -u(xj , sa ""^"^^^ and 

w(a:j_i/2, i") ~ ''^"_i/2- The symmetric spatial difference operator is defined as {Dv")j ~ ''^jVi/2^^j"-i/2 ^^^^ 
(_Du"+^/^)j_i/2 = Uj^^^^ — For the mass and potential term we use a" = (a"^^^^ + a"^"^^^)/2 for 

the integer spacial grid-points and — + '^^'-1/2) ^'-'^ ^^"^ half-integer ones, where a = m,V. 

The space-time stencil is shown in Fig. [T] (a). As one can see, the components u and v are not defined on 
the same space-time grid points but on sub-grids shifted by a half time and space grid spacing. A Taylor 
expansion about the grid points {xj,tn) and (xj_i/2, ^n-i-i/2) shows that the staggered-grid leap-frog scheme 
is second order accurate in space and time. 

The leap-frog time stepping is not self starting. Assume, for example, initial data xp{x,t) given at the 
time-level t". Then the u component has to be propagated from t^ to t^/^. In principle this can be done by 
any time stepping method. We have chosen a symplectic Euler step with half space-time grid spacing, 
which is algorithmically equivalent to the leap-frog scheme, with the exception that all the data (u and v) 
are stored at the same time and space level. Moreover it is only first order accurate (omitting indices for m 
and V): 

^^+^(— ^)^^ + ^^;^-0, ,eZ/2. (22) 
After this first initialization step the w-component is stored only at integer and the w-component at half 



integer spatial points. With a relabeling of the indices, the algorithm Eq. (21 1 can be written as follows: 



with j £ Z, n £ IMq. Here we use the approximations u" « u(xj_i/4, i""-^/''), ~ w(xj.|_i/4, r'+^Z"*) (see 
Fig. [1] (b)). Rearranging of terms leads to the following equations for the explicit recursive update 



z(m- V)At „ 2At/Ax 
2 + i{m-V)At^^ ^ 2 + i{m - V)At 
„+i _ 2 + i{m + V)At „ _ 2At/Aa: 



u-+^ = . 1 



2-i{m + V)At^^ 2-i(m 



(24) 
(25) 



The starting procedure for the leap-frog scheme with a symplectic Euler step of half grid-size is especially 
suitable since symplectic Euler is algorithmically equivalent to the leap-frog staggered-grid scheme and has 
the same dispersion relation for equal ratio r. 

3.2. Von Neumann stability analysis 



For constant coefficients in Eq. (23), Fourier analysis can be used to perform a von Neumann stability 



analysis and to derive the dispersion relation for the whole space problem. A Fourier transform in space of 
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Figure 1: (color online), (a) Staggered-grid scheme for Dirac equation in Pauli-Dirac form with leap-frog time-stepping; (b) 
the algorithmically equivalent (except for starting procedure) symplectic Euler time-stepping. 



Eq. (23) leads to 




. (26) 



Using the definitions ^ = kAx, ji — mAt, v — V At^ and r = At/ Ax the eigenvalues of the amplification 
matrix 

C 2i-u+ti 2ir(e'^-l) \ 

2i+iy-tJ. 2i+i^-n \ 

2r(l-e-'^)(2+»i.-iM) - jJ^^ +iti.i+S.r^ [cos £,-1] I ' > 

are computed by means of 



A± =P/2±^(F/2) -Q , (28) 

where P = tr[G] and Q = det[G], with 

P ^ [2{i'^ - fi^) + 8{1 - + r'^ cos ^)] /N (29) 

and 

Q= [M'-(^-2z)2]/iV, (30) 

where N = fj!^ — {i' + 2i)'^ . A lengthy but straightforward computation yields |A±| — 1 for r < 1. The 
eigenvalues are non-degenerate except for the case fc = /x = and r — 1 (with geometric multiplicity 2), as 
well as the case k = it/ Ax, v = 0, and r ~ 1 (with geometric multiplicity 1). Except for the latter case, the 
scheme is stable under the constraint r < 1, which constitutes its Courant-Friedrichs-Lewy (CFL) condition. 
Below we use a different technique to prove stability which allows one also to identify a functional (related 
to the ?2-iiorm defined on the staggered grid) which is exactly conserved. 



3.3. The energy-momentum dispersion relation 

The Fourier transform of Eq. ( 26 1 with respect to the time variable with a) = LuAt leads to the homoge- 
neous system 

(e'"A B)i^ = , (31) 



with the solutions for r < 1: 



(i± = -^ln[A±] , (32) 
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Figure 2: (color online). Numerical dispersion relation uj(k) for the staggered-grid leap-frog scheme with m = 1 and r = 1 for 
real k (hyperbolically shaped curves) and for imaginary k inside the energy gap (elliptically shaped curve) shown as solid lines. 
The dispersion relation for the continuum problem is shown for comparison using dotted lines: (a) very coarse grid Ax = 1.5 
and k S [tt/Ax, tt/Ax] ; (b) finer grid Aa; = 0.5 where merely an excerpt of the wave number range k £ [7r/(3Aa;), 7r/(3Ax)] is 
shown. 




Figure 3: Dispersion relation u){k) for V = 0, m = 1, and Ax = 0.1, and k e [0, tt/Ax] with r = 1, 0.99, 0.95, 0.9, 0.8, 0.7, 0.6, 0.5 
(from top to bottom). 



where A± is given in Eq. (28). By setting iijV — they reduce to 



hi <^ 1 + ( cos C - l) ± 



r2(cos^ - l) + 1 



- 1 



and with the choice r = 1 they read 

u!± = — i In ^ cos ^ ± y/ cos^ ^ — 1 



i hi ( cos £,±i sin ^ ) = ±^ . 



(33) 



(34) 



Thus for fi^v = 0, and ?' = 1 the hnear dispersion of the continuum problem (Weyl equation) is exactly 
preserved. 



The connection to the phase of the growth factor (i.e. eigenvalues of G) can be established via 

It 1 
ui± = — ln[A±] = — [In |A±| + i arg(A±) + 27rm] = -[arg(A±) + 27rn] . 
r r r 



(35) 



In Fig. [2] the disp ersion relation for a rather coarse grid with Ax = 1.5 is compared to that of a finer one 
using Ax = 0.5, for m = 1, = 0, and r = 1. Clearly, the quality of the numerical dispersion relation 
improves for a finer grid and, for all wave-numbers k € [— tt/Ax, tt/ Ax], it approaches the continuum form 
in the limit Ax and r = 1. In Fig. [3] the dispersion relation for fixed values of ni and V is shown for 



several values of r and Ax = 0.1 . 



3.4- Phase error and gauge invariance 

One can define the phase error for one time-step as 



ephs.sc{k,r,At) = [Wanalytical(fc) ^ ^^nmncnca.l{k, r, At)] At , 



(36) 



where the analytic dispersion relation is ^analytical = ±vr7T? + fc^ and the numerical one is given in Eq. (32). 
Together with Eq. ( 34 ) it follows that Cphaso vanishes in the case of m = ^ = and r — 1 where the scheme 
indeed propagates the solution without error. In Fig. |4] ephaso is shown for different values of m and V on 
a coarse grid At ~ Ax — 0.01 (a) and (b), and on a finer grid At ~ Ax = 10^** (c) and (d). 

Due to gauge invariance, when adding a diagonal term V (constant scalar potential) to the Hamitionan 
Eq. ([T]) one can introduce a new spinor 



Tpit) = '(/;(i)exp ( - iVt) 



(37) 



which fulfills the original equation. Therefore the gauge error per time-step is equivalent to Cphasc due to a 
finite V (see Fig. |4](b) and (d)). 

3.5. Stability analysis within a multiplication technique 

The von Neumann analysis above and the fact that the dispersion relation w(fc) is real for r ^ 1 revealed 
the stability conditions of the scheme for constant mass and potential terms. However, there is another 
technique available which is not based on the Fourier transform and leads to the identification of an "en- 
ergy" functional which is exactly conserved by the scheme, even in the presence of non-constant mass and 
potential terms. The use of this multiplication technique has been inspired by a stability analysis of a 
leap-frog pseudo-spectral scheme for the Schrodinger equation |30j . 

We define the inner product {u,v) := '^jUjVj on P(Z;C) and we will use the notation := {u,u). 
Taking the inner product of Eq. (l20| with (i2"+i/2 + u"-i/2) and taking the real part gives 



1+1/2 



'1/2 



(i:»w",'a"+i/2_^u"-i/2) 



= . 



Eq. (21 1 is multiplied by (w"+^ -|- w") and again the real part is taken to give 



,,"+1 1 



= 



Performing a summation by parts with vanishing boundary terms at infinity gives 



(7^z;",w"+i/2 + u"-i/2) 



-5R 



(£,^t"+l/2^£,^n-l/2^^«) 



Then adding Eq. (|38|) and Eq. (|39f leads to 

-r3fi[(7:>w"+i/2^w"+i) 



,"+1/2 



,ri+l I 



,"-1/2 



{Du 



n-l/2 n. 



(38) 



(39) 



(40) 



(41) 



and one immediately identifies the conserved functional 



,"+1/2 



,."+1 1 



+ r^ 



{Du 



n+l/2 „n+lN 



— const — K 







(42) 
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Figure 4: The dependence of the phase error tphaso the wave vector k S [0, tt/Ax] for y = 0, m = 1, 2, 3, 4, 5 (from bottom 
to top): (a) coarse grid At = Ax = 0.01, (c) fine grid At = Ax = 10"'*; for V = 1,2,3,4,5 (from bottom to top), m = 0: (b) 
coarse grid At = Ax = 0.01, (d) fine grid At = Ax = 10~^. In (d) the functions exceed the plot-range, where the maxima for 
y = (1, 2, 3, 4, 5) are tphasc = (1> 2, 3, 4, 5) * 10"* at fc = 10**7r. 

With this resuh one obtains the stabiUty condition for the scheme by using 



3? 



This gives the estimate 



n+l/2 



< 



,,ri+l I 



2+1/2 



< 



1 - r 



V n , 



(43) 
(44) 



for r < 1. The case r = 1 must be treated separately and one rewrite i?" as 



"+1/2 _ ^n+l/■2^ n+1 

3 )^3 + l/2 



j 



n+l/2 n+1 |2 1 l„ "+1/2 , |2 



"i+l/2l 



Z^l"i+i +V1/21 



Alternatively, when shifting indices, one obtains 



n+l/2 „+l |2 , n+l/2_ „+l ,2 

2^' J" + J-1/2I +2^' ■''+^ J+1/2I 

j j 



(45) 
(46) 

(47) 
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Using i llai + asf < i( Iki + bf + \\ai ~ bf) < I \\ai + bf + ^ \\ai - bf gives 
from Eq. ([46]): u"+^/'^ ^ ^ 

from Eq. (|47]): := XI 



n+l/2 , n+1/2 2 



,,"+1 I ,,"+1 



(48) 



(49) 



This allows us to conclude that the scheme is stable for all r = At/ Ax < 1. Moreover we have identified 
the functional which is conserved by the scheme (see Eq. (42)). In fact, it is conserved for arbitrary time 
and space dependent m,V £ K. 



3.6. Time-reversal invariance 



The time reversal invariance of the scheme can easily be seen in Eqs. (20) and (21). One has to set 
At —At and replace the role of the old and new time-levels n — 1/2 ^ n + l/2 and no n + 1. Then 
one observes that the scheme for the backward propagation has exactly the same form as for the forward 
propagation and concludes the scheme is time-reversal invariant. 



4. Discrete transparent boundary conditions for the staggered-grid leap-frog scheme 

Having discussed the properties of the leap-frog scheme we now turn to the derivation of the associated 
TBCs. Again, we divide the entire space into the computational domain (0,-L), a left semi- infinite exterior 
domain (— oo,0], and a right semi-infinite exterior domain [L,oo). This corresponds to a typical device 
simulation geometry (scattering scenario) in which the nano-device is placed in the computational domain 
and the (macroscopic) contacts are represented by the exterior domains. We make the following simplifying 
assumptions: 

• The initial data ipa^x) = xp((),x) is compactly supported inside the computational domain. 

• In each exterior domain the mass m{x,t) = m and potential l^(a;,t) = V are constant in t and x. 
Both assumption are made for simplicity but can be loosened when needed, as will be discussed below. 



We first Z-transform Eq. (23) in ^-direction on each of the exterior domains, j £ {...,—2,-1,0} and 
j G {J, J+l,J + 2,...}, and solve the resulting finite difference equations explicitly. Using the definition of 



the Z-transform b{z) := Z{bn) = '^b^z its shifting property Z(6„+i) = '^^bn+iz " = '^bnZ 

71—0 n— n— 1 

zb{z) — zbo, and setting 6o = (since the initial spinor is compactly supported on (0,L) ) we obtain 



1 

Ax 





1 

" Ax 
i{m+V) 



Uj{z) 
Vj{z) 

(^+1) 



(50) 



Uj-l{z) 
Vj-liz) 



= . 
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Translation from grid point j to j + 1 therefore is given by 



Ax 1^ 
At 2 

At2 



iV 



i{m+V)Ax 1+z 

\-7? , (m^-\/")(2 + l)" 
At 4 



(51) 



This can also be written as 



where M 



Ax 
At 



(l-z)- 





i(m—V)Ax 



(1 + ^) 



Ax^ 



Ax 1-z I i(rn+V)Ax 1+z 
At z 2 



At2 




a 



fulfills ab = c. Solving the system Eq. ( 52 ) leads to 



Uj+i — (2 + c)uj + Uj-i = 
whose characteristic equation has the roots: 



Tl.2 



C / 

l+2±f +4- 



(52) 



(53) 



(54) 



The same result is obtained for the component v. Since tiT2 — 1, there is one decaying mode (as j oo) 
with |ti| < 1 and one increasing mode with \t2\ > 1. In order to have a solution in ^2(2) one has to choose 
the mode with |ti| < 1. At each boundary only one spinor component couples to the contact region. The 
spinor components u and v, respectively, only need a right and a left TBC (see Eq. ([23| and Fig. [T]) 



voiz) = Ti{z)vi{z) 
uj{z) = ri(z)u,/_i(z) 



(55) 



Then, the inverse Z-transformed boundary conditions are in the form of a convolution in the discrete time 
variable: 

n 
k=0 



(n-fe) fe 



(56) 
(57) 



k=0 



Ti(z) is a non-rational function of z, hence there is no easy way of finding an analytic expression for its 
inverse Z-transformed in general. However, the poles and branch-points in the z-plane, which determine the 
general behavior of their inverse Z-transformed in (discrete) real time n, can be identified. We give a brief 
discussion of these special points for A := At = Ax. 

First, one observes that ti has no poles. For the branch cuts we examine the branch points of the square 
root function - due to zeros of its argument. This leads to the four branch points of ti (Fig. Isl): 



4:(l + iVA) + (m^ -V^)A^ , 2i-(V±m)A 
Zi = -i- , Z2 = -— — ^ . „ and 23,4 = 



4(1 - iVA) + (m2 - F2)A2 



2i+{V±m)A 



(58) 



It can easily be seen that they all lie on the unit circle. These branch points induce damped oscillations 
in time for the convolution coefficients t\^\ The high frequency part (—1)" of these oscillations, which 



arises from zi 



-1, can introduce numerical problems because of subtractive cancellation. It may be 
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advantageous to eliminate this behavior: A multipHcation of ti(z) by (z + l)/z allows one to construct new 
coefScients as a linear combination of the original coefhcients involving the two time steps f^*^) — r'^'^-' -\-t^''~^'> 



(see [211 EH). Then ^ becomes 

n-l 
k=0 

where we used Ti — lim2_>.oo ti{z) = 0. 

For our numerical simulations, the inverse Z-transformation was carried out by performing a power 
series expansion about z = using Mathematica. In Fig. [g] we show the convolution coefficients t^"^ for 
At = Ax = 0.1 and various values of m and V. For the special case m = V = and r = 1 the inverse 
Z-transform can be computed analytically. One gets 



The convolution coefficient then simply is t^""^ = 5", where 5^ is the Kronecker symbol. In the Appendix 
|A]we shall give an explicit, analytic derivation of the convolution coefficients r}""* in the general case. 

Numerical examples 

With the general properties of the scheme established and the associated TBCs identified, we now put 
it to test in challenging numerical applications. 
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.J*- WjJ- '■iiV^ '"'i'i'i- 

50 100 



150 



200 



Figure 6: (color online). Decay of the convolution coefficients r^"' for different values of m and V for Ai = Ax = 0.1: (a) 
m = I and V = 0, (b) m = and V = 1, (c) m = 1 and V = I, (d) m = 1 and V = 2. The darker (blue) color shows 
the real part of the coefficients, the brighter (orange) color shows the imaginary part. The zeroth coefficient r}"' is always 



zero and is not shown in the plots. The first coefficient exceeds the plot-range and has the values: (a) r£ 



(1) 



0.9975, (b) 



.(1) 



0.9925 + 0.0995i, (c) r{^^ = 0.9901 + 0.0990i, (d) t|"' = 0.9682 + 0.1951i 



(1) 



In Fig. [7] wc show the resuhs of a simulation run for a Gaussian wave packet starting out with a mean 
value of fc = 15.92% of k^ax and a standard deviation of ak ~ 1.59% of kmax- A coarse grid Ai = Aa; = 0.05 
is used. The individual figures show the probability density \xjj{x^t)\'^ computed with 



n|2 



i+l/2|2 



"j + i/2l 



3? 



-1/2 



-'j + l 



n-|-l/2N n-fl 



(60) 



which gives the conserved functional Eq. ( |42[ ) when summed over j on an infinite grid. Additionally, the real 
part of the upper component u, and the mass-gap {— energy region between the minimum of the electron 
band and the maximum of the positron/hole band) is shown. One can observe that the wave packet leaves 
the simulation domain without refiections (the numerical error is less than machine precision = 2.2204e-16). 

In Fig. [8] the time dependence of the probability for finding the Dirac particle inside the computational 
interval ||'0(t)|p (:= Eq. (601 summed over j on the computational domain) is shown for this simulation. 

In Fig. [9] we show the reflection at a mass barrier {— spatial region where the energy of the wave packet 
lies inside the mass-gap) with m = 7 and a width of 25 grid-points placed near the center of the simulation 
region. A variation of the mass term can be achieved for Dirac fermions on topological insulators by magnetic 
texturing [34 . The wave packet initially has a mean value of fc = 1.59% of k^ax and a standard deviation 
of (Tfc — 1.99% of kmax, where At = Ax = 0.01. In this tunneling problem it is of profound importance that 
the dispersion relation is correct also for imaginary fc, which is the case for the proposed scheme already for 
relatively coarse grids, as illustrated in Fig. |2] 

Fig. 10 shows a simulation for an initial wave packet with a mean value of fc = 27.06% of k^ax and a 
standard deviation of a^. = 1.99% of kmax moving across a linear potential drop. In this process, the wave 
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X X 

Figure 7: Simulation run for a coarse grid At = Ax = 0.05 and an initial Gaussian wave packet with a mean value of k = 15.92% 
of kmax and a standard deviation of Uf^ — 1.59% of kmax- The figures show the probability density |i/j(x')p (solid lines), the 
real part of the upper component u (dotted lines), and the mass-gap (point-dotted lines). The zero line on the vertical axis is 
shifted by 10, which is the mean energy of the wave packet. With time evolving, the probability HV'II^ for finding the particle 
in the computational domain approaches zero because the wave packet leaves the domain without reflection. 
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Figure 8: Time dependence of the spinor norm of the simulation run shown in Fig. [7] logarithmic scale (loft) and linear 

scale (right). 

number grows beyond the maximum wave number provided by the grid k^ax, where the wave packet is 
under-sampled. The propagation is still correctly described because the dispersion and therefore the group 
velocity is well approximated for all possible wave numbers resolved by the grid. 

4-2. N on- compactly supported initial condition for the spinor and time- dependent exterior potential 
The constraints stated at the beginning of chapter |4] can be loosened: 

• The DTBCs can readily be generalized to the case where the initial wave packet is not compactly 
supported inside the computational domain x £ {0,L). The inhomogeneous boundary conditions can 
be derived by substituting xjj{xr.ht) with xp{xr^i,t) — ■0j„(a;r,;, t) and carrying out the procedure as 
detailed in chapter [ij Here, tp — -i/^j^ should initially have compact support in (0, L). And ■j/'j„ should 
be a solution to the discrete exterior domain problem, e.g. a discrete plain wave (see |31| for details in 
the analogous Schrodinger case). Re-substitution afterwards leads to the result 

n 

n 

< = E H"''^') + -oV • (62) 

fe=0 

The rationale is that the Dirac equation is a linear differential equation which here is approximated 
by a linear finite difference scheme. Therefore one may simply add inhomogeneous boundary terms 
via the superposition principle. 

• The scalar potential V{t) on the exterior domain may depend on time. This is useful when considering 
a time-dependent net potential drop across the system. As for the case of the Schrodinger equation 
a time-dependent exterior potential, e.g. at the right boundary, can be incorporated by treating it in 
the interaction picture which removes V{t) from the Hamiltonian and leads to a gauge (phase change) 
for the spinor Eq. ([g]) |28]. 
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5. Conclusions and future work 

In this paper we have presented a finite difference scheme for the numerical solution of the time-dependent 
Dirac equation in (1+1)D, allowing fully time- and space-dependent mass and potential terms. Owing to a 
combined staggering of the grid both in space and time, unphysical additional Dirac cones are avoided and, 
for the special case of the Weyl equation, the linear dispersion is preserved exactly for all wave numbers 
supported by the grid. In the case of finite mass and/or potential terms the dispersion relation improves 
for all wave numbers, approaching the continuum dispersion relation exactly, when the grid is refined. This 
is a relevant feature when modeling Dirac fermions on a lattice, such as Dirac fermions propagating on 
topological insulator surfaces, graphen, or in quantum spin Hall states. A stability analysis of the scheme 
was performed and a functional, exactly conserved by the scheme, was identified. It provides a valid norm 
for the spinor on the proposed staggered grid. 

Furthermore, we have derived exact DTBCs to close the finite difference scheme for constant mass and 
potential in the boundary regions. With these BCs one is in the position to deal with particle transport 
scenarios and to account for the multi-band nature of the Dirac equation which may cause inter-band 
transfer rather than quantum confinement. For completeness, DTBCs were also derived for the (1+1) Dirac 
(differential) equation in Schrodinger form. 

Using the norm for a measure, numerical simulations of Gaussian wave packets show them leaving the 
computational domain without refiections (with an error below computer precision), thus verifying the 
quality of the DTBCs numerically. 

The assumptions of constant mass and potential on the exterior domain can be loosened. It is however 
desirable that an analytic solution on the discrete exterior domain can be found. The case of purely time- 
dependent exterior potential can be treated by using a proper phase change of the spinor. Initial conditions 
for which the initial wave packet is not compactly supported on the computational domain can be handled 
as well, leading to inhomogeneous terms in the boundary conditions. There are various ways to extend the 
proposed leap-frog scheme to the (2-|-l)D Dirac equation where it retains many of its attractive features 
[32l [33l [35] . The formulation of DTBCs for the (2-|-l)D Dirac equation in the spirit of this paper is the 
subject of future work [35] . 
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Appendix A. Analytic derivation of the convolution coefficients t^"'^ for the general case 



We shall first consider the case m — V — where ( 54 ) reads 



V^z^ - 2fiz + 1 



with fi :— 1 — 2r^. Here, the branch of the square root is chosen such that ti{z) — 0{\z\ ^) as z — >■ oo. 
Using 



Z 



_i j y'z'^ - + 1 



(A.l) 



where P„ denotes the Legendre polynomials (with the convention P_i = P_2 = 0), we obtain 



An) 

'1 - V- ^2^"u ' 2r^"' 2r2 



, n > 0. 



And this simplifies to r}"'' — 5'1 for r = 1 



Next we shall discuss the general case with m, F G R. Here we have 

1 



Ti{z) = 1 + — - — VcVc + 4z, 



where c{z) := zc{z) = az^ + /3z + 7 is a quadratic polynomial given by (52). Using again ( A.l ), both square 
root factors can be inverse Z-transformed. And the explicit formula for the coefficients r^"' would then 
involve a discrete convolution. 

But we shall proceed differently here and rather derive a recursion relation for tI""* (similar as in §3 of 
[36]). A lengthy, but straightforward computation shows that f\{z) :— zti{z) satisfies the inhomogeneous 
differential equation 



c(c + 4z) t[ - [2a^zJ + 3a(/J + 2)z^ + {^^ + 4/3 + 2aj)z + (/3 + 2)7] n = 2z{-f - az"^) . 



(A.2) 



Since all coefficients in ( A.2[ ) are polynomials we shall use the Laurent series of fi, i.e. fi = J2^=o s''^^ z 
withr}"^ A comparison of the coefficients then yields: 

,(0) ^ 1 ,(1) ^ ,(2) ^ /3^+4/3 + 5-a7 



and the exact recursion 

(n+5)a2s("+3) + (2n+7)a(/3+2)s("+2) + (n+2)(/32+4/3+2a7)s("+i)+(2n+l)(/3+2)7s(")+(n-l)72^ = , 
for n > with the convention s^^^^ :— 0. 
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A manuscript showing the extension of the (1+1)D scheme to a (2+l)D scheme which conserves the monotone dispersion 
relation (= having a single Dirac cone) is in preparation. 
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Figure 9: Reflection at a mass barrier with a width of 25 grid-points and a height of m = 7. At = Az = 0.01, the initial 
Gaussian wave packet has a mean value of A; = 1.59% of kmax and a standard deviation of Cfe = 1.99% of kmax- The individual 
figures show the probability density (solid lines), the real part of the upper component u (dotted lines), and the mass 

(dash-dotted lines). The zero for the vertical axis is shifted by 5, which corresponds to the mean energy of the initial wave 
packet. 



21 



iiv(i=i)r= 1 




^ 



-20 L 



50 100 150 200 250 300 350 400 



<> 20 



||V(t=100)f = 1 




50 100 150 200 250 300 350 400 



^ 



;(t=2oo)r= 1 




!> 20 



;(t=3oo)||^= 1 




50 100 150 200 250 300 350 400 



50 100 150 200 250 300 350 400 



Figure 10: Initial Gaussian wave packet with a mean value of fc = 27.06% of kmax and a standart deviation of (7^ = 1.99% 
of ^max meeting a linear potential. A rather coarse grid — Ax — 0.05 is used. The figures show the probability density 
as dashed line, the real part of the upper component u as solid line, and the mass-gap as dash-dotted lines. The zero 
line is shifted by 17 which is the value for the mean energy of the wave packet. 
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